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^ ■ Abstract 
S3 

We propose the variational quantum cavity method to construct a minimal energy subspace of 
wave vectors that are used to obtain some upper bounds for the energy cost of the low-temperature 
excitations. Given a trial wave function we use the cavity method of statistical physics to estimate 

S3' 

the Hamiltonian expectation and to find the optimal variational parameters in the subspace of 



wave vectors orthogonal to the lower-energy wave functions. To this end, we write the overlap 
between two wave functions within the Bethe approximation which allows us to replace the global 
orthogonality constraint with some local constraints on the variational parameters. The method 



is applied to the transverse Ising model and different levels of approximations are compared with 

C3 



the exact numerical solutions for small systems. 
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I. INTRODUCTION 



The variational problem of finding an excited quantum state can be formulated as a 
classical optimization problem constrained by the orthogonality conditions imposed by the 
lower-energy eigenstates of the Hamiltonian. This is a computationally hard task even 
for the ground-state problem in one-dimensional quantum systems. Exact diagonalization 
algorithms are very helpful but they are limited to small systems due to the exponential 
growth of the Hilbert space with the size of system l|. For larger systems one has to resort 
to approximation methods, e.g., the variational quantum Monte Carlo algorithms, to study 
the low-energy states of the Hamiltonian [2 - 4 ] . On the other hand, one can always obtain 
useful insights by studying some exactly solvable mean- field models j^]. 

The cavity method of statistical physics provides the framework and machinery we need 



to deal with random optimization and constraint satisfaction problems js-9|. In particu- 
lar, the method provides a local message-passing algorithm to study single instances of a 
problem on random and finite-connectivity graphs within the hierarchy of replica symmetry 



breaking (RSB) approximations jsl, [lo|, |Uj. For instance, the cavity method can be used 
to minimize the Bethe estimation of computationally difficult cost functions like the Hamil- 
tonian expectation in a variational quantum problem; first the expectation is computed in 
the replica symmetry (RS) approximation to obtain an energy function of the variational 
parameters. Then, the parameters are considered as statistical variables to study their sta- 
tistical properties within a higher level RS approximation, with equations that resemble 
the one-step RSB ones. These equations are exploited in a message-passing algorithm to 
optimize over the variational parameters. There are other applications of the cavity method 
to quantum systems that work with the density matrices or the path integral representation 
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of the quantum partition function 



In this paper we use the variational quantum cavity method [17| to find approximate so- 
lutions for the exited-states of the transverse Ising model. The main point of this study is to 
write the overlap between two wave functions in the Bethe approximation, which is asymp- 
totically exact as long as the trial wave functions can be represented by classical systems of 
locally tree-like interaction graphs. This allows us to replace the global orthogonality con- 
dition with some local constraints on the variational parameters and the cavity messages in 
the Bethe expression for the overlap. In summary, given an appropriate trial wave function, 



we evaluate the Hamiltonian expectation in the subspace of wave functions orthogonal to 
the lower-energy states of the Hamiltonian and find an estimation of the optimal variational 
parameters, all within the Bethe approximation and implemented by a local message-passing 
algorithm. 

Note that, usually, for large and disordered systems we do not have the exact ground- 
state of the system. The best that we can do is to find a good variational wave function that 
represent the low-energy state of the system. Obviously, an excited state that is orthogonal 
to an approximate ground state does not necessarily provide an upper bound for the excited- 
state energy. There are, of course, other ways to find the excited states that do not rely on 
the orthogonality to a priori known ground state, but we found the above procedure more 
amenable to the cavity method that we are going to use in this paper. Nevertheless, any 
orthogonal set of quantum states defines a subspace of wave vectors that according to the 
Courant min-max theorem [18] can be used to find some upper bounds for the Hamiltonian 
eigenvalues. 

In the next section we give the definitions and write the general equations. Section II III 
is devoted to the mean-field approximation, where we take product states to represent the 
low-energy eigenstates of the Hamiltonian. Then, in Sec. |IV]we consider more complicated 
trial wave functions including the two-body Jastrow interactions. The numerical results are 
presented in Sec. [V] and the concluding remarks are given in Sec. EH 



II. DEFINITIONS 

Consider the transverse field Ising model with Hamiltonian H = — Yl(ij)es q Jij a i a j ~ 
Y2ihi a i with index i = 1,...,N that labels the sites in the quantum interaction graph 
S q . The a^' y,z are the standard Pauli matrices. In the following we will work in the a z 
representation with orthonormal basis \a) = |eri<72 • - • ctn}- Given a trial wave function 
|\P(P)) = So-VKffi £)\^) depending on a set of variational parameters P, we write the 
Hamiltonian expectation as 

(*(P)\H\^(P))=J2m^P)\ 2 E(a), E{a) = £ eij {a h a s ) + £ ei (a), (1) 
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where e^o*, oj) = —JijOiOj and 



^(-)--^|x:^§(-vri-)|. (2) 

Depending on the trial wave function we obtain different expressions for the non-diagonal 
terms e;(a). 

We consider fi(a;P) = \^{^P)\ 2 as a probability measure in a classical system and 
compute the above average quantities within the Bethe approximation. For P_) oc 
Hcie^c ( l ) a((I da "i Pa), the classical measure is given by P) oc Y[ a \4>a((L da ", Pa)\ 2 with the set 
of classical interactions S c = {0 a (a Sa )|a = 1,... , A}. Here 9a is the subset of variables that 
appear in <ft a . Similarly we define di as the subset of interactions that depend on <7j. The 
cavity marginal /ij_ >a (a'j) gives the probability of having <jj in absence of interaction term 
4> a . Similarly we define ^^(cj) as the probability of having Cj in absence of the other 
interactions involving <jj. The equations governing these cavity marginals are: 

fJ>i-+a(<Ti) OC I j Hb^CTi), (3) 

OC IM^Pall 2 ]J Hj^afa)- (4) 

a aa \' j£da\i 



These are the belief propagation (BP) equations 19] that can be solved by iteration starting 
from random initial cavity marginals. In the replica-symmetric approximation we assume 
;here is a fixed point to the BP equations describing the single Gibbs state of the system 



7] 



. The above cavity marginals are enough to obtain the Bethe estimation of the average 
energy. 

Let us denote by \^n(P_ n )) the nth trial wave function that minimizes the average energy 
(^n(E n )\H\^ n (P n )) conditioned on the orthogonality constraints (^ n (P n )\^m(P m )) = 
for m = 0, ...,n — 1. The corresponding classical systems are represented by measures 
/x n (a;P n ). In the following we are going to satisfy the orthogonality constraints within the 
Bethe approximation, 

(y n (P n )\y m (P m )) = ^<(ffi£ n )Va<i;£ m ) - e-^^^'-Sw^) = o, (5) 

where AF{, AF a , and AF ia are the free energy changes by adding variable node i, interaction 
node a, and link (ia) to the complex measure v n , m &P n ,P m ) oc ^nfc P n )i J m{<I.', ET) oc 



n a 4>* a (c[_ , Pa)<t>a{<L°"i PT)- These quantities are given by 



0, 



e 



e n ^(^)' ( 6 ) 

E O0a(£ 5a ; PD II ( ? ) 

y]vjr+g{(Ti)VarM(<7i), (8) 



where the cavity marginals z/j_ i>a (crj) and z/ a _j,j(<jj) satisfy the BP equations for the complex 
measure v n ,m(<I.', P" - , P m ) ■ Thus, to have orthogonality it is enough to have e~ AFa = for 
some a. This defines a constraint on the parameter P™ given P m . 

In summary, to estimate the average energy and to satisfy the orthogonality constraints 
we need to know the BP marginals of the classical measure fi n and {z/ njm |m = 0, . . . ,n — 
1}. In addition, we have to choose the set of constrained parameters A n = {a m \m = 
. . . , n — 1} for the orthogonality constraints e~ AFam = 0. In this paper we will use a greedy 
strategy to construct A n , by choosing the parameters that at least locally minimize the 
energy expectation. Finally, the problem of minimizing over the variational parameters is 
considered classical statistical physics problem 



Z « = E E 1bp II w^-^K (9) 

E" n n ,{iyn, m } m=0,...,n-l 



where for f5 opt — > oo the Gibbs measure is concentrated on the optimal parameters. The 
indicator functions Ibp and I„ >m ensure that the messages fi n , {v n ,m} satisfy the BP equations 
and the states n and m are orthogonal. Starting from n = 0, one can find the other states 
one by one after solving the above optimization problem. 

One may find an approximate solution to the above problem by a two-stage algorithm: 
Given {P m |m = 0, . . . ,n — 1} and an arbitrary set of the constrained parameters A n , we 
run BP to find the set of marginals {v n ,m,\ m — 0, . . . , n — 1}. These are used to fix the 
constrained parameters in A n to satisfy the orthogonality constraints. Then we minimize 
the average energy (E(a))^ n over the remaining parameters. The above two stages are 
repeated to converge the algorithm. 



III. THE MEAN-FIELD APPROXIMATION 



Let us start with the mean-field (MF) approximation, where the trial wave functions are 
represented by the product states: 

ip&B) oc JJe BlCT % (10) 

i 

with complex parameters Bi. This results to ej(<7j) = —hie~ 2B ? ai cos(2B-(Ji) and the follow- 
ing classical measure /i(cr;.B) oc n i e lB * i ' Ti - By superscripts R and I we mean the real and 
imaginary part of the parameters. Given the above measure we find 

e^i „ = -h — ' 11 
cosa(2B^) 

(eijiat, <Tj))p = -Jij tanh(2Sf ) tanh(2Bf ). (12) 



We see that for non-negative hi the average energy is minimized by setting B\ = 0. 
Therefore, as long as we are interested in the ground state, we can set the imaginary parts 
to zero. In this case, the variational problem reads 

B 

The cavity marginals of the parameters in the Bethe approximation are 

M^iBi) oc e-A**^" J] lj2 e ~ Me ' k) " M ^( B k))- (14) 

i.-^h j \ B, J 



For (3op t oo and scaling M^^Sj) = g-^^^C- 8 ;) we g nc j ^ e m i nsum equations 



M ir +j(B i ) = (e i ) li + m&{(eik)r + Mk-*i(B k )}. (15) 



B k 
k£di\j 



The equations are solved by iteration starting from random initial messages. After each 
iteration we subtract a constant from the messages to have mine. Ai i-^j(Bi) = 0. Then we 
find the optimal parameters by minimizing the local minsum weights, 



B? = argmin \ (e^ + Y] min {<e ii ) A , + M^B,)} j> ■ (16) 



Note that in the MF approximation the orthogonality condition reads 

(^ n {B n )\^m(B m )) oc n cosn ( 5 r + B?) = 0, (17) 
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thus, it is enough to have cosh(S™* + Bf 1 ) = for some i. This means Bf R = —Bf 11 * and 
Bf 1 — Bf = 7r/2. Therefore, the nth excited state can be obtained by the following set of 
constraints: 



ln,m = { m m (18) 

Bf - Bf = f , 



for m = 0,...,n — 1. The index z m is chosen to minimize the local energy (ej m ((Tj))^. Then, 
one can use the same minsum equations as above to minimize over the remaining parameters. 
In this way we can find at most N orthogonal product states of minimum energies E n . 

Given the states \^ m {B m )) for m = 0,...,n, we can easily compute the Hamiltonian 
matrix elements H mm i = (^ m (B m )\H\ i if m '(BJ n )), which for arbitrary parameters read 



it _ TT ( cosher + BY ) 

mm ' 11 V (cosh(2Ef i? ) cosh(2Bf R ))V 2 



X 



£ Ja tanh(^ + Bf) tanh(^ r + Bf) - ^ fr ^gC %\ \ ■ (19) 



One can diagonalize the Hamiltonian in the subspace spanned by the above states to obtain 
the eigenvalues A m . Then, using the min-max principle, we know that A ra is an upper bound 
for the nth eigenvalue of the Hamiltonian. Indeed, we found that in this case A n ~ E n 
for large N as the off-diagonal matrix elements H mm > decay exponentially with the size of 
system. 

Notice that instead of imposing the orthogonality exactly we could ask for an exponen- 
tially small overlap (^ n {B n ) \^ m (B m )) < e N by demanding 



cosher + sr) <f£ij 



{cosh(2Bf R ) cosh(25 i mfl )) 1 /2 



for all i, which can easily be imposed in the above minsum equations. Indeed to apply the 



min-max theorem we do not need a set of orthogonal states 18]; according to the theorem, 
the nth eigenvalue is given by mins n+1 max|^) 6 g n+i: ^|^) = i(^|f/"|^i) where 5 n+1 is any subspace 
of dimension n + 1. 
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IV. BEYOND THE MEAN-FIELD APPROXIMATION 



We can do better than the MF approximation by adding the local two-body or Jastrow 
interactions [20] to the trial wave functions: 

4>ij(ai,aj), (21) 

with 

&(<7i) = e B ^\ fa^cTi, a,) = e K ^. (22) 
For simplicity we are going to assume £ c = £ q . As a result, we obtain 



-h t e- 2B ^~^^ ' 1K " n n - cos(2fl/<7i + 2K ijWi)- ( 23 ) 

j£di 



The average energy is computed with respect to the following classical measure 

Hi I0i(°"i)| 2 Yl(ij)ee c l<M°"i>°"j)| 2 - To estimate the average energies (ey^, er,-))^ and 
(ej(<jj, a*)) M we need the following local marginals 

liij{(T h (Tj) = ^-|0 ii (ai,a i )|Vi->i(^)/"j^i(° r j) ) (24) 



^i,di{vuv dl ) = -J—lM^l 2 II l0ii(c r i,o-j)lVj^i( (J j)> (25) 

^i,di 

given in terms of the cavity marginals 



fii-+j((7i) oc |0j(cTj)| 2 J [ | ^ ^^(a^cTfc)! 2 ^^^^) I • (26) 

k&di\j \ (Tk / 

For Bj = we can simplify the equations by taking the symmetric (or paramagnetic) 
solution of the BP equations. This is of course exact for a tree classical interaction graph 
£ c . Then, the average local energies are given by 



j€di 



eiM, aj))^ = -Jij tanh(2AT£). (28) 



The resulting minsum equations are 



AW#ii) = trr ™ I (e,) M + « e «>^ + >W^)) } , (29) 
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and the optimal couplings are estimated by 

K% = arg min {(e^ + M^K^) + M^K^)} . (30) 
The orthogonality condition for two symmetric states reads 

(^> n {K n )\m rn {K m )) oc eZ^c {K Z* +K ™ )aiaj oc Y[ cosh(^* + K%) = 0, (31) 

— {ij)e£c 

using the symmetric solution of the BP equations. Thus, to have orthogonality we need 
cosh(i^™* + K™) = for some link (ij). That is, for the nth excited state we have the 
following set of constraints: 

f K? R : + K™* =0, 

U-n,m — \ yoZ,) 

K — K ■ = — 

for m = 0, ... , n — 1. The link (i m j m ) is chosen to minimize the local energy (ei m j m ) Mn . 
Then, one can use the same minsum equations as above to minimize over the remaining 
parameters. The number of orthogonal states that we can find in this way is limited by 
the number of the coupling parameters in the classical system which for a tree structure is 
N — 1. The Hamiltonian matrix elements are given by 

cosh(fT™* + K$) 



H mm i 




cosh(2K" f lR ) cosh(2ir" f l '- R )) 1 /2 



X 



.S^*r + <)-E».jn[(^^jj- (33) 

Similarly, we can make more general orthogonal states \^ n (B^, IC 1 )} and |\& n (B_ m , K™)) 
by choosing Bf such that for the complex measure v nm oc e^^ Bi +Bi "> ai+ ^{ij)ee c ( K ij + K i j) <T ^j 
we have e AFi = for some i. This gives 



,2(BJ»+BJ") _ _ 



( e 



The node « can be chosen in a greedy way to minimize (ei(<7j, <7 9l )) Mn in the mean-field 
approximation. Here it is more difficult to minimize the average energy {E(a)) fln , which 
depends not only on the variational parameters but also on the BP cavity marginals fj Ji ^j(a i ). 
More precisely, the minsum equations read 

Mi v i7\ <; .// (; ) min <{e i ) fl + (fak)? + M k ^i(K ik , fj, ik )) > , (35) 
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FIG. 1. The ground-state and excited-state energies (Eq and E{) for the transverse Ising model 
on a random tree of size TV = 20 with random Gaussian couplings Jjj of mean zero and variance 
one in a uniform transverse field hi = h. The exact results are compared with the upper bounds 
that are obtained by a minimal energy subspace spanned by product states (MFoMF), symmetric 
states (SoS), and tree states (ToT). 

where /i^ = (fj,i->j, fJ>j-n)- Note that the minimum in the right hand side is conditioned on 
satisfying the local BP equations. The reader can find more details in Ref. 



V. NUMERICAL RESULTS 



Let us start with a small system to compare the above approximations with the exact 
results for the ground and excited states. We take a random tree with random Gaussian 
couplings Jij of mean zero and variance one in uniform transverse fields hi = h. This system 
displays a phase transition from the ordered phase for h < h c with a nonzero Edwards- 
Anderson order parameter q = 'Yl li (o'l) 2 /N to a disordered phase for h > h c . 

As Fig. [1] shows, we obtain better ground-state energies with the product and symmetric 
states in the ordered and disordered phases, respectively. While the product states allow for 
a nonzero magnetization, the symmetric states have by definition zero magnetization and 
therefore more appropriate to represent the disordered ground state. Indeed, an estimate 
of the transition point can be obtained by comparing the ground-state energies that are 
computed by the product and symmetric states. Figure [2] displays these energies for the 
one-dimensional transverse Ising model with ferromagnetic couplings. Nevertheless, we see 
in Fig. [1] that for small system sizes we always obtain better energies for the excited state 
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FIG. 2. The ground-state energy density (Eq/N) for the one-dimensional (ID) transverse Ising 
model with uniform and ferromagnetic couplings (Jjj = 1) obtained by the product states (MF) 
and the symmetric states (S) in the thermodynamic limit. 

by a minimal energy subspace of the product states. In Fig. [1] we also display the results 
obtained by more general trial wave functions having a tree structure defined by S c = £ q . 
As expected, we obtain much better upper bounds by introducing both the variational 
parameters B{ and K{j. Figure |3] shows the results for a larger number of spins. Here, in the 
disordered phase we find better upper bounds for the excited-state energy by the symmetric 
trial wave function. This is due to the presence of very small couplings that reduce the cost 
of the orthogonality constraint. 

Finally, we present the MF results for the transverse Ising model in a two-dimensional 
square lattice. In Fig. H]we compare the upper bounds for the ground and excited states 
with the exact ones in a small lattice. 

VI. CONCLUSION 

We generalized the variational quantum cavity method to study low-temperature excita- 
tions of quantum systems within the Bethe approximation. We constructed orthogonal sets 
of minimal energy quantum states, where the Hamiltonian matrix elements can be computed 
exactly to obtain some upper bounds for the Hamiltonian eigenvalues. For more general trial 
wave functions we have only an approximate estimation of the Hamiltonian matrix elements 
but the estimation is expected to be asymptotically exact as long as the trial wave functions 
are represented by locally tree-like classical interaction graphs. And finally, the method can 
also be extended to include some appropriate global interactions in the trial wave functions 
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FIG. 3. The ground-state and excited-state energy densities {Eq/N and E\/N) for the transverse 
Ising model on a random tree of size N = 1000 with random Gaussian couplings Jij of mean zero 
and variance one in a uniform transverse field hi = h. We compare the upper bounds that are 
obtained by a minimal energy subspace spanned by product states (MFoMF), symmetric states 
(SoS), and tree states (ToT). 




FIG. 4. The ground-state and excited-state energy densities {Eq/N and E\/N) for the transverse 
Ising model on a two-dimensional (2D) square lattice of size N = 4 x 5 with random Gaussian 
couplings of mean zero and variance one in a uniform transverse field hi = h. The exact results 
are compared with the upper bounds that are obtained by a minimal energy subspace spanned by 
product states (MFoMF). 



that are essential to deal with the fermion sign problem 211 ] . 
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